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High spectral and spatial resolution, mid-infrared fine structure line observa- 
tions toward two ultracompact HII (UCHII) regions (G29.96 -0.02 and Mon R2) 
allow us to study the structure and kinematics of cometary UCHII regions. In 
our earlier study of Mon R2, we showed that highly organized mass motions ac- 
counted for most of the velocity structure in that UCHII region. In this work, 
we show that the kinematics in both Mon R2 and G29.96 are consistent with 
motion along an approximately paraboloidal shell. We model the velocity struc- 
ture seen in our mapping data and test the stellar wind bow shock model for 
such paraboloidal like flows. The observations and the simulation indicate that 
the ram pressures of the stellar wind and ambient interstellar medium cause the 
accumulated mass in the bow shock to flow along the surface of the shock. A 
relaxation code reproduces the mass flow's velocity structure as derived by the 
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analytical solution. It further predicts that the pressure gradient along the flow 
can accelerate ionized gas to a speed higher than that of the moving star. In the 
original bow shock model, the star speed relative to the ambient medium was 
considered to be the exit speed of ionized gas in the shell. 

Subject headings: interstellar medium - ultracompact H II regions - G29.96 -0.02 
- Mon R2 IRS1 



1. Introduction 

Ultracompact HII (UCHII) regions form when massive OB stars ignite inside dense 
molecular clouds and high energy UV photons ionize the surrounding neutral material. They 
have small sizes (<0.1 pc), high electron densities (>10 4 cm -3 ), and high emission measures 
[EM = Jnin e dl > 10 7 pccm~ 6 ) (Wood & Churchwell 1989b). The ionization increases the 
temperature (~10 4 K) and number density (H 2 — > 2p + + 2e~) causing the gas pressure inside 
the regions to be at least two orders of magnitude higher than in the surrounding molecular 
gas. UCHII regions should then expand at approximately the speed of sound in the ionized 
gas (~ 10 kms -1 ) until they reach pressure equilibrium with the surrounding material or 
break out of their parent molecular clouds. Because of this expansion, the density and 
emission measure both should drop with time. Once the emission measure drops below 
10 7 pccm -6 , these HII regions will no longer be considered ultracompact. 

A longstanding puzzle regarding UCHII regions involves their numbers and their life- 
times. If HII regions expand at 10 km s _1 , they should remain ultracompact for only ~ 10 4 yr. 
Since this is <1% of the lifetime of an OB star, the number of UCHII regions should be <1% 
the number of OB stars. Radio interferometry observations confirm that the ratio of IRAS 
far-infrared (FIR) flux densities provides an efficient way to select embedded OB stars from 
the IRAS data (Wood & Churchwell 1989b; Churchwell 1990; Garay et al. 1993; Kurtz et al. 
1994; Miralles et al. 1994). However, counts of IRAS FIR sources with the spectral character- 
istics of UCHII regions and the optically visible O stars in the solar neighborhood show that 
OB stars spend ~ 10 — 20% of their main-sequence lifetime embedded in molecular clouds 
(Wood & Churchwell 1989a). Apparently, the lifetimes of UCHII regions are an order of 
magnitude larger than predicted by the classical pressure-driven spherical expansion model. 

The variety of observed UCHII morphologies, including spherical, shell-like, cometary, 
and irregular, also needs to be explained. An understanding of the kinematics of the ionized 
gas will help provide an explanation of both the morphologies and the lifetime of UCHII 
regions. It has been suggested that long lifetimes could result from some kind of containment 
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or gas replenishment mechanism (Hollenbach et al. 1994; Dyson et al. 1995; Redman et al. 
1996; Williams et al. 1996; Redman et al. 1998). Various mechanisms would have different 
effects on the morphologies and the kinematics of UCHII region, suggesting that observations 
could distinguish between these mechanisms. 

The well-organized overall appearance of cometary UCHII regions make them good 
objects in which to study UCHII region kinematics and morphology. One way to explain 
these objects is with a bow shock model, which suggests the cometary structure of some 
UCHII regions is the result of the supersonic motion of OB stars with high speed winds 
through molecular clouds (Hughes & Viner 1976; Wood & Churchwell 1989b; van Buren 
et al. 1990; Mac Low et al. 1991; van Buren & Mac Low 1992). Swept-up ambient gas and 
stellar wind material accumulate where the ram pressures associated with these two mass 
flows balance, resulting in a shell-like structure and a surface flow. This model may not 
explain the lifetime of all UCHII regions, but could provide a part of the solution to this 
problem. 

Because ultracompact HII regions are formed inside dense molecular clouds, the ex- 
tinction toward these objects usually is very high, with a typical line of sight extinction of 
A v = 30 — 50 (Hanson et al. 2002). Thus, optical techniques used in diffuse HII region 
observations are inappropriate for studies of UCHII regions. Low extinction at infrared and 
radio wavelengths makes these wavelength bands more suitable for studies of UCHII regions. 
Hydrogen radio recombination line observations toward UCHII regions have been carried 
out extensively (Garay et al. 1985; Kim & Koo 2001; Araya et al. 2002). These observations 
reveal that thermal motions cannot account for the linewidths of many UCHII regions. How- 
ever, significant thermal broadening makes it hard to study the velocity structure in these 
regions using hydrogen recombination lines because of hydrogen's low atomic mass (Jaffe & 
Martin-Pintado 1999; Sewilo et al. 2004; De Pree et al. 2004, and references therein). The 
thermal line width of hydrogen recombination lines is Avfwhm— 21-4 kms -1 for gas with a 
temperature of 10 4 K. This large thermal linewidth means that heavier ions are better probes 
of bulk motion; the thermal broadening is only about 4.8 kms -1 for Ne + ions at the same 
temperature. Mid-infrared ionic fine-structure lines have been used to probe the structure 
and excitation of HII regions (Beck et al. 1981; Lacy et al. 1982; Takahashi et al. 2000), but 
generally with too low spectral resolution to study gas motions. To study kinematics and 
distinguish organized motions from turbulence and thermal broadening of heavy ions, one 
needs ~ 5 kms -1 velocity resolution. 

We have begun a program of high spectral resolution observations of mid-infrared fine 
structure line emission from UCHII regions. By mapping a small sample of UCHII regions 
of different morphologies at high spatial resolution in [Ne II] A12.8 /im, [Ar III] A9.0 jum, 
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[S III] A18.7 yum and [S IV]A10.5 /im, we can study the kinematics and physical conditions 
in UCHII regions in order to better understand massive star formation and the relationship 
between these stars and their surrounding environments. Among these lines, the [Ne II] line 
is particularly bright and is present in HII regions with a broad range of excitation, so it 
is well suited for kinematic analysis. We presented high spectral resolution observations of 
one cometary UCHII region, Mon R2 IRS1 (Mon R2 hereafter), in Jaffe et al. (2003). In the 
current paper, we present [Ne II] observations of another cometary UCHII region, G29.96 - 
0.02 (G29.96 hereafter). We also examine bow shock models in some detail and compare 
them to the observations of both G29.96 and Mon R2. 

In section 2, we describe our mapping and data reduction methods. We present our 
[Ne II] line emission observations toward G29.96 in section 3. In section 4, we describe the 
kinematics of the bow shock model and our relaxation method for simulating the formation 
of the surface flow. We also introduce an additional acceleration, due to the pressure gradient 
in the ionized gas along the flow. In section bf 5, we compare the model predictions with 
the [Ne II] observations of G29.96 and Mon R2. Finally, we discuss the existing problems of 
models for cometary UCHII regions in section 6. 



2. Observations 
2.1. Instrument 

The [Ne II] observations of G29.96 were carried out with TEXES (the Texas Echelon 
Cross Echelle Spectrograph, Lacy et al. 2002) on the NASA 3 meter Infrared Telescope 
Facility (IRTF) on Mauna Kea, in June, 2001. TEXES is a high resolution (R < 100,000) 
spectrograph operating at mid-infrared (5-25 /zm) wavelengths. The slit, with a 1.4" width 
and a 11.5" length, was oriented north-south on the sky. With this slit, we achieved a spectral 
resolution of ~4 kms -1 or R = 75,000 at 12.8 fim. Each pixel along the slit was about 0.36" 
on the sky, which is a little smaller than half of the diffraction limit (0.88") of the IRTF 
at 12.8 fxm. To map the regions of interest we stepped the telescope west to east across 
the objects without chopping. Multiple, partially overlapping scans were needed to cover 
the entire region. For G29.96, each scan had a length of 20" with a step size of 0.4". For 
Mon R2, each scan was 45" long with 0.7" steps (Jaffe et al. 2003). Nine overlapping scans 
were made on G29.96, and 20 scans were made on Mon R2. Because the scans overlapped, 
the integration time at each point in the map is approximately 18 seconds for G29.96 and 8 
seconds for Mon R2. 
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2.2. Data Reduction 

The spectral-spatial datacubes produced from the scans were first reduced with a cus- 
tom Fortran reduction program that performs general purpose corrections, such as the 
correction of optical distortions, flat fielding, bad pixel masking and cosmic ray spike re- 
moval (Lacy et al. 2002). Radiometric calibration was obtained from measurements of 
an ambient temperature blackbody before each set of scans, giving intensities in units of 
ergcm~ 2 s _1 sr _1 (cm -1 ) -1 . The uncertainties in the intensity are mostly systematic, and are 
probably ~ ±20%. The wavelength calibration was obtained from sky emission lines, and is 
accurate to ~1 kms -1 . IDL scripts were used to do the sky background subtraction, multiple 
scan cross-correlation and combining, and datacube manipulation. Positions safely off the 
object at both ends of each scan allow us to interpolate the sky emission at each step posi- 
tion. In doing the sky background subtraction, we assumed that the sky background varies 
linearly with time during the course of a scan, based on the brightness of pixels off the ob- 
ject. After sky emission subtraction, multiple scans of the same region were cross-correlated, 
shifted, and added to make a complete datacube for the object. 



3. Results 

We show the integrated [Ne II] line emission map for G29.96 in Fig. 1. The overall 
cometary shape is apparent and is symmetric about an axis at a position angle ~ 70° east of 
north. The brightest emission forms an arc perpendicular to the symmetry axis. To the west 
of the arc, the line emission drops rapidly from the peak value to the background level within 
4". On the east side of the bright arc, fainter emission extends over 13" to the east-northeast. 
The change of the emission level is much more gradual and less uniform. The east edge of 
the region is much fuzzier and broken by a faint "lane" into two parts. The southern portion 
extends the curve of the emission arc out ~ 6" and ends at a clump of ionized gas of size 
~ 3" x 3". From our velocity channel maps (Fig. 2), we can see that it is the continuous 
extension of the arc. 

Velocity channel maps of G29.96 (Fig. 2) show additional details. A cavity is present 
at medium and lower velocity channels, from ~ 82 km s" 1 to ~ 94 km s -1 . The fainter 
structures on the east side of the region and the brighter arc form an almost complete ring 
around the cavity, although the symmetry axis of the elliptical ring is at a somewhat larger 
position angle (~ 120°) than that of the bright arc (~ 70°). A similar ring is seen in channel 
maps of Mon R2 (Jaffe et al. 2003), except that the major axis of the ring in Mon R2 is nearly 
perpendicular to the symmetry axis of its bright arc. A compact emission peak centered at 
Vlsr ~ 108 kms" 1 can be seen ~ 5" northeast of the peak of the arc in G29.96. The [Ne II] 
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line is generally broad, reaching ~ 40 kms" 1 near the bright arc. The emission is typically 
bright near the velocity of the ambient cloud (Vlsr ~ 98 kms -1 , obtained from single dish 
observations of rotational transitions of CS molecules, Mart in- Hernandez et al. 2003, and 
references therein), but the line is significantly blueshifted ahead of the emission peak. The 
line center reaches a Vlsr ~ 82 km s -1 at the west edge of the region 

In Fig. 3, we plot [Ne II] line contours on the top of the VLA 2 cm continuum image of 
G29.96 (Fey et al. 1995). We first aligned the two maps, assuming they are emitted by the 
same ionized gas. We then convolved the radio map with our beam. The close resemblance 
of the maps after convolution suggests that the much sharper edges seen in the VLA image 
are the result of its higher spatial resolution (< 0.56"). Comparing a cross-cut along the 
symmetry axis (Fig. 4) confirms that [Ne II] and radio free-free flux are both proportional 
to the emission measure of the region. Radio recombination line observations at 0.62" and 
4 km s -1 resolution showed similar arc-like structure (Wood & Churchwell 1991). They also 
showed that the H76a line is broader in front of the radio continuum emission arc, where 
the velocity gradient is also the highest. 

Martin-Hernandez et al. (2003) obtained a narrow-band H 2 1-0 S(l) filter map and K 
band (2.07 — 2.19/iin) spectra of the region with a resolution of R=8,000 along a 120" long 
slit along the symmetry axis. Their map shows many local near-infrared emission clumps 
together with the proposed ionizing star of the region lying ~ 2". 3 northeast {p. a. = 64°) 
of the emission peak. We mark the star position with an asterisk in our map according to 
this offset (Fig. 1). Their Bvy linewidth ranges between ~42 kms -1 and ~62 kms -1 and 
shows significant blueshift in front of the emission arc, which is consistent with our [Ne II] 
observations. The fact that Br^y line is systematically broader than [Ne II] along the slit 
supports the idea that [Ne II] is a better tracer of the ionized gas motions. 

4. Bow Shock and Relaxation - A Physical Model for Parabolic Flow 

In this section, we present an extension of the analytic solution for bow shocks using a 
numerical relaxation technique. We then use the resulting gas flow as a template with which 
we can compare the observed kinematics of G29.96 and Mon R2. 

4.1. Analytical Solution 

Bow shock models have been explored to explain the formation of cometary UCHII 
regions by several authors (van Buren et al. 1990; Mac Low et al. 1991; van Buren & Mac 
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Low 1992; Wilkin 1996). They describe a situation where a star with a strong stellar wind 
and high UV luminosity moves inside a dense molecular cloud. In the frame of reference 
of the star, the stellar wind and ambient material collide and create a stationary shock 
region approximately paraboloidal in shape in front of the star. The standoff distance r is 
determined by pressure balance at the apex of the shell: 

p w vl + n w kT w = p a vl + n a kT a , (1) 

where p w = is the stellar wind density at the distance r and p a = n a p is the density 

of the ambient medium. The corresponding number densities are n w and n a . Where the 
stellar wind and ambient medium collide near head-on, the ram pressures are much greater 
than the gas thermal pressure, thus the standoff distance can be expressed as: 



/ M w v w , 

by neglecting gas thermal pressure (Wilkin 1996). The gas pressure becomes non-negligible 
far downstream from the star where the shell surface is roughly parallel to the stellar motion, 
resulting in reduced ram pressure from the ambient medium. 

For physical and mathematical simplicity, most models assume that the stellar wind and 
ambient medium material are well-coupled inside the shock region, and radiative cooling is 
very efficient on the time scale of the ionized gas recombination time. The resulting shell 
is "momentum supported", rather than "energy supported", and relatively thin compared 
to its scale. Mass and momentum are conserved and transported within the shell. This 
momentum conserving assumption provides the possibility of investigating the kinematics 
analytically. Wilkin (1996) derived the formula for the shell's shape by applying mass and 
momentum conservation and neglecting the gas thermal pressure on both sides of the shell, 
as well as the tangential acceleration caused by pressure gradients within the shell: 



r(0) = t o esc 9 sj 3(1 — 9 cot 9) , (3) 

where r(0) = r is the standoff distance, and 9 is the angle between the point on the shell 
and the apex. He also derived expressions for the surface density and the tangential velocity 
of the shell. 



4.2. Relaxation 

The analytical method has no numerical uncertainties, but necessarily leaves out some 
details. Once we include more terms in the momentum equation, it can no longer be in- 
tegrated analytically. Assuming a bow shock is a steady state configuration, we can apply 
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a relaxation method to the mathematical description of the system. Using the relaxation 
method allows us to include more physical processes than the analytical method of Wilkin 
(1996), while being easier to formulate and less time-consuming to calculate than a full 
hydrodynamic approach. 

We assumed a cylindrically symmetric bow shock shaped as an approximately paraboloidal 
shell. A grid with fixed angular size is created on the surface of the shell and each cell has 
the same angular size along both the azimuthal and the polar directions, 86 = 5<p = 0.5°. 
Due to the symmetry of such a geometry, we only need to calculate quantities in one strip 
of cells along the polar direction. We assume that the gas in a cell moves as a single fluid. 
That is, we assume rapid cooling behind the shock where the ambient and stellar wind ma- 
terial mix and form a uniform shell, and we neglect shear motion between front and back 
sides of a cell. Mass flow, gas velocity, and position calculated from the analytical formula 
(Wilkin 1996) are assigned to each cell as initial conditions (which were shown not to affect 
the results). Then, we let the program iteratively adjust these parameters by applying mass 
and momentum conservation. 

During relaxation, the physical and geometric properties of each cell are recalculated 
to replace the previous values in each step until the whole system achieves convergence. In 
each step, the mass flow from one cell (i) into a neighboring cell (i+1) is given by: 

Sm l = 5m i ~ 1 + 5ml + bm\ (4) 

where 5m l ~ l is the mass flowing out of the previous cell. The mass flows from the stellar 
wind and the ambient medium are 5m % w and Sm l a . They are given by: 

Sml = j^M w t (5) 
5m\ = -p a (v a .S*)r (6) 

Sfl is the solid angle subtended by the cell element, M w is the mass loss rate of the stellar 
wind, S l is the outward pointing normal vector with magnitude equal the surface area of the 
cell i, v a is the velocity vector of the ambient material, and r is the time step of the iteration. 
Under the same condition, we also have the formula for total momentum flow from the cell: 

5 P * = Sp*' 1 + 5m^ w + Smivi + (P* - P*)SV + $ P (7) 

where P* and P l w are the gas pressures of the external media. The orientation of the normal is 
perpendicular to the direction of tangential velocity, V 1 = 5p' l /5m\ $ P is the extra pressure 
term resulting from the tangential pressure gradient inside the shell. We include it in the 
relaxation calculation after testing our model for the case without pressure acceleration. The 
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distance of the cell from the star is calculated from the law of sines: 

r* +1 = Sin( ^ } r' (8) 
sin(/3* -69) ' { } 

where 58 is the angular separation between cells and (3 is the angle formed by the radius 
and the velocity vector. Because we assumed that the system will eventually reach a steady 
state, we did not try to solve a non-steady problem. Our calculation simply proceeds and 
finds a converging solution. 

The shape of the shell and the tangential velocity along the shell, calculated by the 
relaxation method neglecting the contribution of the gas pressure gradient, are shown (dashed 
line) in Fig. 5. For this illustration, we set the star's speed equal to 20 km s -1 . For 
comparison, the analytical solution (Wilkin 1996) is also shown (dotted line). We also 
show an overall shape of the calculated shell in Fig. 6. In this later calculation, we include 
gas thermal pressure contributions from both sides of the shell by assuming a temperature 
of 50K for the ambient medium and 10 4 K for the free flowing stellar wind, which gives 
a more accurate but only slightly different solution for the tail portion of the shell than 
the analytical solution. We can see that the discrepancy between our iterated result and 
the analytical solution (Wilkin 1996) is smaller before we add pressure acceleration to the 
relaxation solution. Calculations show that, in the rest frame of the star, gas starts to move 
with zero tangential speed from the apex of the paraboloidal like shell and accelerates until 
it reaches the star's travelling speed, 20 kms -1 , at the end of the shell. This acceleration is 
due to the accumulation of mass and momentum, which are provided by swept-up ambient 
material and the stellar wind. In later sections, we will show that this simple picture may not 
be totally correct since the effect of the pressure gradient within the shell can be significant. 

Fig. 5 also shows the surface density and the particle number density along the shell. The 
surface density increases, almost linearly, from its value at the apex of the shell. The particle 
number density, derived from the shape of the shell and the normal pressure components, 
drops rapidly from its value at the apex. It drops by a factor of four by the position where 
the shell passes the star. In steady state, gas in the shell should be in appoximate pressure 
equilibrium with the average of the normal components of the pressures on the two sides. 
The difference between the two normal components causes a pressure gradient across the 
shell which causes its centripetal accelaration. From the pressure in the ionized gas, we can 
calculate the density: n = Pram -^ p " as ; where we take the ionized shell's temperature as 10 4 K. 
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4.3. Pressure Gradient Acceleration 

Most existing solutions, both analytical (Wilkin 1996) and numerical (Mac Low et al. 
1991), neglect the effect of the gas pressure gradient along the compressed shell. The ex- 
ception is Comeron & Kaper (1998), who make a numerical hydrodynamic calculation, but 
concentrate on the case of a runaway O star moving through the diffuse ISM. In addition 
to the centripetal acceleration caused by the unbalanced external pressures, the variation in 
the pressure along the shell results in a pressure gradient that accelerates the gas along the 
shell. The momentum deposited in a cell by this effect is given by: 

pi—l pi+1 

$ P = AV = yP^r (9) 

Where P l and \/P l are the pressure and the pressure gradient in the ith cell. A 1 and 
V 1 are the cross section and the volume of the cell. We derive the thickness of each cell 
from the emission measure, using the ionization-recombination equilibrium equation, and 
the density, calculated from the balance between gas pressure and ram pressure assuming 
stellar parameters appropriate for G29.96 (section 4.4.1). We then assume that the thickness 
of the given cell is constant. If the central star is unable to ionize the whole shell, the thickness 
of the shell is slightly greater than that of the ionized layer, with the neutral gas in a thin, 
dense region just outside of the HII region. The thickness increases from the apex of the shell 
to the tail region because the density drops along the shell. The results of the relaxation 
are shown with a solid line in Fig. 5. The analytical solution is shown as a dotted line for 
comparison. 

From the plots in Fig. 5, we can see that the change to the shape of the shell caused 
by including the pressure gradient is very small. Thus the change in the ionized gas number 
density is also small, because it is calculated from the ram pressure normal component, which 
directly relates to the shape of the shell. In the region close to the apex, this change can 
be neglected. A larger change is seen in the tangential velocity plot. For the parameters 
assumed, which includes a partially ionized shell, the velocity near the head of the bow shock 
increases so that the gas speed is ~ 1 km s _1 higher than without the pressure gradient 
acceleration at 9 ~ 1 radian. Without the pressure gradient, the maximum velocity that 
gas can achieve is the speed of the star. With the additional acceleration provided by the 
pressure gradient, the gas in the shell can reach higher velocities. The tangential speed 
finally drops as the shell picks up more mass from the ambient medium and the pressure 
gradient acceleration effect decreases toward the end of the shell due to the drop of stellar 
wind ram pressure. With the parameters we have chosen, the maximum speed that gas can 
achieve under the pressure gradient acceleration is 1.3 km s -1 higher than the speed of the 
star. Up to this point, we have assumed that the ionized gas and the swept-up neutral gas 
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are coupled, so they share the momentum carried by the stellar wind and ambient medium. 
When including the effect of a pressure gradient, we assume that the force resulting from 
this effect is also shared by ionized and neutral gases. In fact, most of this force should be 
exerted on the ionized gas only, since a pressure gradient causes a force per unit volume and 
the ionized gas fills most of the volume. If the ionized gas can slip past the neutral layer, 
it should reach a higher velocity than what we present in this work. The resulting speed 
should be similar to that in a fully ionized shell situation. 

If the optical depth of the shell is small for ionizing radiation, the ionizing photons will 
ionize the whole shell or penetrate the shell and ionize the ambient medium beyond it. In 
this "ionizing photon leaking" situation, the pressure gradient acceleration effect will reach 
its maximum because the force is proportional to the thickness of the ionized shell. Our 
calculations show that the gas speed in the shell can reach 1.25 times the star speed in the 
range 10 - 20 km s _1 , and there is no drop in speed for angles larger than 2.4 radians seen 
in Fig. 5(b). 

In this work, we only show the results of the single layer bow shock model. Including 
multiple layers to our model is a logical next step. Furthermore, we leave out the pressure 
gradient in the hot post-shock gas on the inner side of the shell from our model. Since it 
is so hot (~10 7 K), it should cool very slowly. It is even harder to include the contribution 
from this part of HII regions in the model. 

The pressure gradient accelerates gas inside the shell, thus decreasing the surface density. 
It does not affect our flux maps because the brightness of each cell in our model depends only 
on the number of ionizing photons, as long as the shell is only partially ionized. However, 
the change in the surface density will affect the thickness of the neutral part of the shell, 
thus affecting the appearance of corresponding molecular flux maps. A kinematic study of 
the neutral component associated with ionized bow shock structures might prove useful. 

If pressure gradient acceleration is ignored, the spatial and kinematic structure of a 
modeled bow shock depend on two parameters: the speed of the star and the stand-off 
distance. The dependence of the stand-off distance on the stellar mass-loss rate, stellar wind 
terminal speed, stellar speed, and the ambient medium density is given in Eq. 2. As is 
discussed above, inclusion of the effect of the pressure gradient in the ionized gas has little 
effect on the spatial distribution, but changes the velocities much like an increase of the 
stellar speed would. A set of models for UCHII regions in G29.96 and Mon R2 is given in 
Table 1, where plausible values, based on other observations (see references in the table), are 
given for unconstrained parameters. 
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5. Comparison with Observations 

Our relaxation models of bow shocks provide the shape of the shell and the Doppler 
shift at each point on the shell. We calculated the [Ne II] emission at each point on the 
shell, assuming ionization-recombination balance in the shell, taking into account the angle 
at which the ionizing radiation hits the shell. We ignore any dust that might be present in 
the region. 

The model ionized shell was tilted relative to the line of sight to improve the agreement 
with the data, sampled and remapped onto the sky with a spatial grid 0.2" x 0.2", and then 
convolved with TEXES's velocity resolution of 4 km s _1 and a model of TEXES's beam 
of half-maximum radius r ~ 0.8". For G29.96, we convolved with an additional turbulent 
linewidth of 10.0 km s -1 to improve the agreement with the observations. In Mon R2, we 
convolved the data with the the narrowest line width ~8.8 km s -1 (Jaffe et al. 2003). Our 
observing sampling grid is much coarser than our modeling grids in the region corresponding 
to our observations, so the error caused by the quantization of the modeling grids is small. 

5.1. G29.96 -0.02 

Four model parameters are constrained with varying uncertainties by our observations: 
the standoff distance (ro), the tilt of the symmetry axis from our line-of-sight, the speed of 
the star relative to the ambient medium and LSR velocity of the star (or equivalently of the 
molecular cloud). In addition, the stellar ionizing luminosity and spectral type, and the neon 
abundance combine to determine the [Ne II] brightness. Other parameters, notably those 
from which the standoff distance can be calculated, can be determined from other observa- 
tions of our objects or estimated from typical O-star properties. Parameters determined in 
this way are given in Table 1. 

We rotate the observed flux datacube 26° (assuming p.a.=64°, Martin-Hernandez et al. 
2003) counter-clockwise in order to orient the axis of the shell horizontally in the figures. We 
also position the observed flux map so that the position of the ionizing star matches that in 
our model. A cross at (0, 0) in the map (Fig. 1) indicates the position of the ionizing star. 
The spatial distribution of the [Ne II] emission is reproduced better if the motion of the star 
in the model is within 45° of the plane of the sky, but a bigger angle tends to give a better 
fit in position-velocity diagrams. We chose to tilt the shell 50° away from us in our model. 

As discussed in section 4.3, we assume that the ionized gas and the neutral gas in 
the shell are moving together, and we pick parameters so that the shell is only partially 
ionized. A stellar speed of v* ~ 20kms~ 1 is picked to fit the observed range of velocities. 
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An LSR velocity of the star of V^^sr = 104 kms x , or V am b,LSR = 89 kms 1 is needed 
to fit the velocity offset observed. However, observations of the nearby molecular material 
(Churchwell et al. 1990; Cesaroni et al. 1992; Afflerbach et al. 1994; Olmi & Cesaroni 1999; 
Lumsden & Hoare 1999; Martin-Hernandez et al. 2003) give Vam^^sR = 92 — 100 kms -1 , 
with most numbers near 98 kms -1 . A stellar speed of v* ~ 10 kms -1 would allow better 
agreement between the model velocity offset and the molecular observations, but would not 
fit the observed [Ne II] velocity range unless the pressure gradient acceleration is larger than 
in our model. 

The simulated [Ne II] flux map and P-V diagrams on cuts through various positions 
in the model are shown with corresponding diagrams from observed data in Fig. 7-8s. The 
accelerating, paraboloidal like flow produced in the bow shock model successfully matches 
many global features in the observed data. The model produces the limb brightening at 
the head region, although the fit would be improved by using a smaller tilt away from the 
observer. The curvature of the bright ridge is also fit well. Cuts perpendicular to the shell axis 
show similar central velocity shifts, spatial and spectral ranges, and overall shape, including 
a ">", in both the observations and the model (Fig. 7). From cut to cut, the curvature of the 
">" changes because the motion of the shell with the star gives more redshift to gas closer 
to the apex than to down-stream gas. Our line-of-sight passes through the ionized shell 
twice for positions near the symmetry axis. Given the orientation of the UCHII region, the 
far side is closer to the apex, contains higher density gas, and dominates the line emission. 
The near side is farther from the apex and has lower density gas. At the edges of the cuts, 
we also see more lower density gas. The ">" forms because gas closer to the apex is less 
blue-shifted, compared with gas farther from the apex. Because the cuts are perpendicular 
to the symmetry axis, the model P-V diagrams are necessarily symmetric. Asymmetries 
in the observations are also small. In cuts parallel to the symmetry axis (Fig. 8), the P-V 
diagrams show a "7" like pattern, with a curved leg and a rather flat top "arm" , indicating 
a large velocity gradient at the head and less velocity change towards the tail. The line is 
also broader at the head region. These are predicted by the model too. In our velocity plot 
(Fig. 5), we can see that more than 80% of the change in the tangential velocity happens 
in the first ir/2 radians. The line is broader because the scale length for velocity change is 
smaller in the head region. The remarkable similarity of the models and the observations is 
a clear indication for a large-scale paraboloidal like flow in this UCHII region. 
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5.2. Mon R2 IRS1 

Our observations of the compact HII region in Mon R2 IRS1 are presented in a previous 
paper (Jaffe et al. 2003). With a 24" diameter shell and a bright southeast ridge, [Ne II] 
emission in Mon R2 shows complex and broad velocity structure. We speculated that we 
were looking at a kinematic pattern in which material flows from the bottom to the rim of 
a bowl-like feature. Here, we try using the bow shock model to interpret the observations, 
since the object also has a cometary appearance. In order to show the similarities between 
the kinematics of a shell-like flow structure and the observed data, we choose parameters so 
that the pressure gradient acceleration is negligible. The results of two models with different 
standoff distances are shown in Figs. 9-14. In the first model, we try to match the ionizing 
star position to that in near-infrared observations that are good to < 1" (Yao et al. 1997). In 
the second model (Table 1), we use a bigger standoff distance, which gives a better fit to our 
[Ne II] observations. In both cases, we tilt the shell 20° toward the observer in the simulation 
for a good fit to the observed P-V diagrams. As in the G29.96 -0.02 case, we derive a V arrihtLS R 
different from the value we found in the literature for Mon R2. An additional ~ 8 km s _1 
redshift is needed to shift the center of the line to the value in the observational data. It is 
interesting to note that Mon R2 requires a redshift and is tilted toward us while G29.96 is 
tilted away and requires additional blueshift. This may indicate the pressure gradient may 
have a bigger effect on the gas acceleration than we calculate and that a more sophisticated 
model might fit the data with a smaller stellar velocity. 

Qualitatively, the model agrees with the data. The observed morphology appears 
cometary, with a bright arc and a fainter tail, although the tail ends more abruptly than in 
the model. The P-V diagrams show two peaks over most of the region, as is expected where 
our line of sight passes through the front and back sides of a shell. 

If the predominant motion were due to expansion of the shell, rather than flow along 
the surface, the velocity splitting would increase only gradually moving from the edge of the 
shell into the center. As in G29.96, the motions in Mon R2 are consistent with gas flow along 
a shell. However, there are a few facts which make Mon R2 HII region hard to fit into the 
bow shock model. First, there is a compact broad-lined region near or just inside of the apex 
of the observed shell, which is not predicted by the model. It is most apparent in Fig. 11 
and Fig. 14, where it produces the central ridge in the P-V diagrams. This component is 
also shown in our previous paper (Jaffe et al. 2003, Fig. 7). This source may be a result of 
the shock front overtaking a dense clump in the molecular cloud, or, it could result from an 
instability in the front, as is seen in the hydrodynamic models of Comeron & Kaper (1998). 
The second difficulty is that a larger standoff distance is required to make the curvature of 
the shell the same as in the data. In the first cut of Figs. 9 and 12, where our observed 
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P-V diagrams only show one component, our model P-V diagrams show two components. 
In addition, the shell appears to be closed on the back side. The outermost contours on the 
images are nearly circular, although the rim is much brighter at the bottom of the map than 
at the top. The model predicts that the shell is always open toward the tail. 

In our figures 9,10 for both the model and the data, cuts parallel to and near the shell 
axis show two components. At the ridge of the shell-like region, these two components are 
connected by a broad line with a width up to 40 kms -1 . The relative strengths and the 
separation of the two components changes with the position of cuts. Normally, the blue- 
shifted component is seen farther toward the tail because the tilt of the shell puts the denser 
part of the near side of the shell farther from the head than the far side. The single broad 
line in the P-V diagram of the cut made in the center of the ridge is probably caused by the 
additional source there. 

The observed P-V diagrams of the cuts perpendicular to the shell axis (Fig. 11,14) 
demonstrate the gradual change of a circular structure, which is typical for cuts across a 
rotationally symmetric shell. If we neglect the broad-line component at the peak of the flux 
map, our model and observed P-V diagrams agree well. Compared to the observations, the 
model with a smaller standoff distance predicts a smaller spatial span (Fig. 11). Once we 
increase the standoff distance (Fig. 14), we no longer have this problem. A density gradient 
along the symmetry axis of the shell might be able to explain the curvature of the shell, 
although there is at present no concrete evidence for such a gradient. 

6. Discussion 

We have shown that the overall velocity structure behind the cometary shape of UCHII 
regions in G29.96 -0.02 and Mon R2 is formed when ionized gas flows along a paraboloidal 
like surface. Using the bow-shock model, we can reproduce this kind of structure. We found 
that the model provides a good qualitative explanation of our observations of G29.96 and 
Mon R2. Both sources have morphologies similar to that predicted by the model, at least 
near the apex of the shell. The observed line profiles, as seen in the P-V diagrams, are 
generally double-peaked, indicating that the gas is swept up into a thin shell, so each line of 
sight passes through two surfaces with different Doppler shifts. The two Doppler components 
reach a maximum separation near the shell apex, indicating that the gas accelerates along 
the shell, and that the dominant motion is along the shell rather than radial. 

There are quantitative differences between the model and the data, especially in the 
case of Mon R2, as discussed above, indicating that a complete description of these sources 
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will have to be more complicated than our simple model. Nevertheless, we think that the 
bow shock model must be essentially valid. Now we consider whether this model resolves 
any of the problems associated with the study of UCHII regions, and how the model might 
be improved. 

The current bow shock model assumes that the stellar wind material and the ambient 
medium mix well and cool efficiently after they collide, so their mass and momentum con- 
tribute to the swept-up layer and the velocity of the ionized gas is that of the single swept-up 
layer. However, this momentum- conserving thin shell assumption is not well justified. Be- 
hind the shock, stellar wind material is collisionally heated and ionized. This gas will have 
a temperature > 10 7 K and does not cool efficiently, keeping the layer thick. This high tem- 
perature layer forms between the stellar wind and the photonionized gas layer and makes the 
mass and momentum exchange between the two layers difficult. The thin shell assumption 
is invalid in this situation. The photonionized layer is also evaporated through conduction 
with the extremly hot shocked stellar wind gas, which makes predicting the kinematics more 
difficult. There are many instabilities which can inflate the shell to make the thin shell 
assumption invalid. Possible candidates include the transverse acceleration instability (TA1) 
(Dgani et al. 1996) caused by the acceleration of the flow normal to the surface of the shell 
and the non-linear thin shell instability (NTSI) (Vishniac 1994; Hueckstaedt 2003) caused 
directly by the collision of isothermal flows. These instabilities will disturb the shell and 
create sub-structures with scales comparable to the thickness of the bow shock. Under such 
conditions, the momentum conservation assumption will not hold. 

Among the other models proposed to explain cometary UCH11 regions, the champagne 
model is the most interesting. The model applies when a massive star is found in a region 
with a large density gradient, such as the edge of a dense molecular cloud. The resulting H1I 
region will expand supersonically away from the high-density region, or simply break out of 
the edge of the molecular cloud and cause the ionized gas to stream out of the opening in 
response to a large pressure gradient. The classic champagne model without a stellar wind 
(Tenorio-Tagle 1979; Bodenheimer et al. 1979; Yorke et al. 1983) can explain the cometary 
shape of Hll regions but, because the ionized gas fills up the bubble and the main pressure 
gradient is along the density gradient, the champagne model has difficulty accounting for 
limb brightening and for the line profiles observed. Adding a stellar wind to the champagne 
model will probably produce limb brightening (Comeron 1997). Without ram pressure of 
the external gas, the champagne model with a stellar wind tends to produce a bigger shell 
in the clouds with the same density and temperature as the clouds we investigated with the 
bow shock model. 

The main differences between the bow shock and the champagne flow are seen in the 
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kinematic properties of the ionized gas (Garay & Lizano 1999). First, the bow shock model 
predicts that the velocity gradient is steeper in the head than in the tail, which can be seen 
easily in our P-V diagrams, whereas the champagne model expects the largest gradient in the 
tail, where the fractional pressure gradient also reaches the highest value (Comeron 1997). 
Second, the champagne model predicts that the line widths are broader in the tail region 
because the gas has higher speed there and the gas motion is less parallel, while the bow 
shock predicts broader line widths near the apex because the gas gains more momentum 
there and moves cylindrical-symmetrically along the surface (in the frame of reference of the 
star). Finally, ionized gas near the apex of the cometary structure is at rest with respect to 
the ambient molecular gas in the champagne model, whereas it moves with the star and/or 
the shock front in the bow shock model. In all of these respects, our observations of G29.96 
and Mon R2 are in better agreement with the bow shock model. 

A central problem about UCHII regions is their lifetimes. The number ratio of UCHII 
regions and OB stars requires that UCHII regions should have an average lifetime ~ 20% 
of that of OB stars (Wood & Churchwell 1989a). Any model of UCHII regions has to be 
able to account for this. We have observed over a dozen of UCHII regions. They all have 
large velocity range, usually over ~ 20 km s _1 . In the bow shock model, the velocity range 
of the line is directly connected to the star speed relative to the ambient medium. If 20% of 
UCHII regions have a bow shock-like structure, as shown in Wood & Churchwell (1989a,b) 
and Kurtz et al. (1994), the same percentage of OB stars should move supersonically through 
molecular clouds. Although evolved OB stars can be accelerated up to 200km s -1 through 
the association ejection or supernova explosions (Blaauw 1993) and form bow shock-like 
structures around them (van Buren et al. 1995; Noriega-Crespo et al. 1997; Kaper et al. 1997), 
high speed OB stars are rare in OB associations. The velocity dispersion of OB association 
is generally small, only a few km s" 1 (Jones & Walker 1988; Tian et al. 1996). Even if the 
molecular gas initially made the potential a bit deeper, stellar speeds of 15 — 20 km s _1 
are extremely improbable as part of a normal distribution. Including pressure gradient 
acceleration helps to reduce this requirement. Our calculations suggest that the needed 
speed is still too high to solve the problem. This disagreement indicates that the stellar 
speeds may be less than in our model. Due to the simplicity, our model could underestimate 
the acceleration of the ionized gas along the shock front. A full hydrodynamical treatment 
would help determine the true effect of a pressure gradient. 

We are grateful to Alan Fey and Ed Churchwell for letting us use their VLA 2 cm 
data (Fey et al. 1995; Wood & Churchwell 1989b). We thank Gregory Shields for help on 
the ionization model. We also need to thank NASA IRTF staff for their help on observa- 
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Fig. 1. — Integrated line flux map of [Ne II] line toward G29.96. The (0,0) postion marks 
the flux peak of the observation, corresponding to 18M6 m 03 s .92, -02°39'21".9 ( J2000). The 
asterisk marks the position of the ionizing star from Martin-Hernandez et al. (2003). Con- 
tours are drawn at 70%, 50%, 35%, 25%, 17.5%, 12.5%, 9% and 6% of the peak value 
(1977 ergs cm -2 s -1 sr -1 (cm -1 ) -1 ). 
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Fig. 2.- 



Channel maps of G29.96. Corresponding Vlsr is indicated in the plots. 
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Fig. 3. — Contours of [Ne II] line map of G29.96 (solid line, peak value = 
1977 erg cm -2 s _1 sr _1 (cm -1 ) -1 ) overplotted on VLA 2 cm continuum observation (dotted 
line, peak value = 0.078 Jy/beam with a 0.56" x 0.49" beam and a 100MHz bandwidth Fey 
et al. 1995). The VLA map has been smoothed with the ([Ne II]) beam profile of 1.6" 
FWHM (peak value = 0.037 Jy/beam after the smoothing). Contours are plotted at 5%, 
10%, 20%, 30%, 40%, 60%, 80%, 95% of the corresponding peak values. Coordinates are 
the offsets from the emission peak. The straight dashed line (p.a.=64°)indicates where flux 
is taken for the flux distribution plot in Fig. 4. 




RA Offset (arcsec) 

Fig. 4. — Fluxes along the symmetry axis (shown in Fig. 3 with the dashed line) of G29.96 - 
0.02. The dotted line indicates the smoothed radio 2 cm flux, while the solid line indicates 
flux of [Ne II]. The curves are normalized to have equal integrals. 
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Fig. 5. — Model for G29.96 -0.02. (a) The shape of the shock, plotted against the angle 
from the apex (6). (b) The tangential velocity V t along the surface, (c) The shape of the 
shock. X is the symmetry axis of the shell. Offsets are measured from the position of the 
ionizing star. The length unit is the standoff distance (ro). (d) The tangential velocity Vt 
along the shock surface, (e) Surface density (a). (f)The particle number density N along the 
bow shock. The results from Wilkin (1996) are shown with dotted lines. The model without 
the pressure gradient is shown with dashed lines. The model with the pressure gradient is 
shown with solid lines. Notice that the results of the model without the pressure gradient is 
extremly close to the results of Wilkin (1996). 
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Fig. 6. — The calculated shape of the bow shock without the pressure gradient contribution 
and the corresponding tangential velocity vectors in G29.96 -0.02 using parameters of model 1 
in Table 1. The cross at (0, 0) indicates the position of the star. r , the standoff distance 
from the star to the apex of the shell, is the length unit in the plot. The dotted line indicates 
that the line of sight from the lower left to the upper right, form 40° with the symmetry axis 
of the shell. 
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G29.96 [Nell] P-V Diagram Model [Ne II] Model P-V Diagram 
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Fig. 7. — Position- Velocity (P-V) diagram of G29.96 -0.02 [Ne II] observation and that of 
model 1 (Table 1). Contours in the flux map are drawn at 50%, 25%, 12.5%, 6.25%, 3.12%, 
1.56% and 0.78% of the peak value. Contours in the P-V diagrams are drawn in the similar 
way, but with contour levels separated by a/2. A cross at (0, 0) indicates the location of the 
star according to the model. Same for the following plots. 
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Fig. 8. — Position-velocity diagram of G29.96 -0.02 [Ne II] observation (top panel) and 
model 1 (bottom panel). 
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MonR2 [Nell] P-V Diagram Model P-V Diagram 
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Fig. 9. — UCHII region in Monoceros R2 and model 2 with the standoff distance matching 
that in Yao et al. (1997) (Table 1). [Ne II] map and the position-velocity diagrams along the 
corresponding cutting lines. Two asterisks and one cross mark the position of IRS1 (0, 0) 
and IRS2 (0, 12.4) from Yao et al. (1997) and the position of the ionizing star (0, 0) used in 
the model. 
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MonR? [Nell] P-V Diagram Model P-V Diagram 
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Fig. 10. — Same as in Fig. 9. 
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Fig. 11. — Same as in Fig. 9. 
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P-V Diagram 
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Fig. 12. — HII region in Monoceros R2 and model 3 with a longer standoff distance than 
that in Yao et al. (1997) (Table 1). [Ne II] map and the position- velocity diagrams along the 
corresponding cutting lines. Two asterisks and one cross mark the position of IRS1 (0, -2.45) 
and IRS2 (0, 14.8) from Yao et al. (1997) and the position of an ionizing star (0, 0) according 
to the model. 
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MonR? [Nell] P-V Diagram Model P-V Diagram 
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Fig. 13. — Same as in Fig. 12. 
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Fig. 14. — Same as in Fig. 12. 
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Type 


Distance 


Radial 


Medium 


Mass loss 


Wind 


Travelling 


Distance 


The Sky 


Angle 








Velocity 


Density 


Rate 


Speed 


speed 












D 


Vls.r 


n a 


M w 


v w 


v a 


r 










(kpc) 


(kms- 1 ) 


(10 s cm" 3 ) 


(10- 6 M Q yr- 1 ) 


(kms" 1 ) 


(kms" 1 ) 


(10 16 cm) 






G29.96 -0.02 


06 a 


6 a 


98 b 


0.5 c 










2.27" d 




Model 1 


05.5 


6 


89 


0.11 


8.0 


1000 


20 


19.6 


2.2" 


50° 


Mon R2 IRS1 


09 ° 


0.95 f 


10.3 g 


0.3 h 










4.35" i 




Model 2 


BO 


0.95 


10.3 


0.244 


3.6 


1000 


27 


6.61 


4.4" 


-20° 


Model 3 


BO 


0.95 


10.3 


0.238 


8.0 


1000 


27 


9.95 


6.8" 


-20° 



Note. — We estimated the spectral type of the central star for our models from our [Nc II] line luminosity, assuming optically thin in the whole 
region. Free parameters of our models are the cloud LSR velocity (Vlsr), the ambient medium density (n a ), the stellar mass loss rate (M w ), the 
stellar wind speed (v w ), the stellar speed relative to the ambient medium (v a ) and the tilt angle. For models, the size on the sky is the distance 
from the emission peak to the ionizing star after the tilting and the beam convolution. Positive tilt angles represent tilting the symmetry axis 
away from the observers. The parameters are chosen so that the pressure gradient acceleration effect is negligible. 

a Pratap ct al. (1999) 

b Afflerbach et al. (1994) 

c Morisset et al. (2002) 

d Martm-Hernandcz ct al. (2003) 

c Downes ct al. (1975) 

f Racine & van den Bergh (1970) 

SMeyers-Rice & Lada (1991) 

h Montalban et al. (1990) 

'Yao ct al. (1997) 



